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ABSTRACT 

We use a set of cosmological simulations combined with radiative transfer calculations 
to investigate the distribution of neutral hydrogen in the post-reionization Universe. 
We assess the contributions from the metagalactic ionizing background, collisional 
ionization and diffuse recombination radiation to the total ionization rate at redshifts 
z = — 5. We find that the densities above which hydrogen self-shielding becomes im- 
portant are consistent with analytic calculations and previous works. However, because 
of the intensity of diffuse recombination radiation, which peaks at the same density, 
the transition between highly ionized and self-shielded regions is smoother than what 
is usually assumed. We provide fitting functions to the simulated photoionization rate 
as a function of density and show that post-processing simulations with the fitted 
rates yields results that are in excellent agreement with the original radiative transfer 
calculations. The predicted neutral hydrogen column density distributions agree very 
well with the observations. In particular, the simulations reproduce the remarkable 
lack of evolution in the column density distribution of Lyman limit and weak damped 
Lya systems below z = 3. The evolution of the low column density end is affected 
by the increasing importance of collisional ionization with decreasing redshift. On the 
other hand, the simulations predict the abundance of strong damped Lya systems to 
broadly track the cosmic star formation rate density. 

Key words: radiative transfer - methods: numerical - galaxies: evolution - galaxies: 
formation - galaxies: high-redshift - intergalactic medium 



1 INTRODUCTION 

A significant fraction of the interstellar medium (ISM) in 
galaxies consists of atomic hydrogen. This makes studying 
the distribution of neutral hydrogen (HI) and its evolution 
crucial for our understanding of various aspects of star for- 
mation. In the local universe, the HI content of galaxies is 
measured through 21-cm observations, but at higher red- 
shifts this will not be possible until the advent of signifi- 
cantly more powerful telescopes such as the Square Kilome- 
ter ArrajQ However, at z < 6, i.e., after reionization, the 
neutral gas can already be probed through the absorption 
signatures imprinted by the intervening HI systems on the 
spectra of bright background sources, such as quasars. 

The early observational constraints on the HI column 
density distribution function (HI CCDF hereafter), from 
quasar absorption spectroscopy at z < 3, were well de- 
scribed by a single power-law in the range Am ~ 10 13 — 
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10 cm" " 2 (|Tvtlerlll987l b Thanks to a significant increase 
in the number of observed quasars and improved obser- 
vational techniques, more recent studies have extended 
these observations to both lower and h igher HI column 
densities and to higher redshifts (e.g. , iKim et alj 2002: 
iPeroux et~al"1l2005l; lo Meara et al.ll2007l; iNoterdaeme et all 
20091: iProchaska et all 120091: IProchaska fe Wolfe] 120091 : 



O'Meara et al.ll2012l ; lNoterdaeme et al]|2012h . These studies 
have revealed a much more complex shape which has been 
described using several different power-law functions (e.g., 
IProchaska et alj|20ld : lO'Meara et al.ll2012r i. 

The shape of the HI CDDF is determined by both the 
distribution and ionization state of hydrogen. Consequently, 
determining the distribution function of HI column densi- 
ties requires not only accurate modeling of the cosmologi- 
cal distribution of gas, but also radiative transfer (RT) of 
ionizing photons. The complexity of the RT calculation de- 
pends on the HI column density. At low HI column densi- 
ties (i.e., Ahi < 10 17 cm -2 , corresponding to the so-called 
Lyman-Q forest), hydrogen is highly ionized by the meta- 
galactic ultra-violet background radiation (hereafter UVB) 
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and largely transparent to the ionizing radiation. For these 
systems, the HI column densities can therefore be accurately 
computed in the optically thin limit. At higher HI column 
densities (i.e., Am > 10 17 cm -2 , corresponding to the so- 
called Lyman Limit and Damped Lyman-a systems), the 
gas becomes optically thick and self-shielded. As a result, 
the accurate computation of the HI column densities in these 
systems requires precise RT simulations. On the other hand, 
at the highest HI column densities where the gas is fully self- 
shielded and the recombination rate is high, non-local RT 
effects are not very important and the gas remains largely 
neutral. At these column densities, the hydrogen ionization 
rate may, ho wever, be strongly affected by the loc al sources 
of ionization (jMiralda-Escudel 120051 ; [Schavc 2006 ; Rahmati 
et al. i n prep.) . In addition, other processes like H2 for- 
mation i|Schavdl200l"bl ; iKrumholz et al.ll2009bl ; lAltav et all 
120111 ) or mechanical fee dback from young stars and / or 
AGNs l|Erkal et alj|2012i ). can also affect the highest HI col- 
umn densities. 

Despite the importance of RT effects, most of the pre- 
vious theoretical works on the HI column density distribu- 
tion did not attempt to model RT effects in detail (e.g 
Katz et all 1 19961; iGardner et~ai]|l997t lHaehnelt et all 1998 



„ .Gardner et al. Ill997l; Hac hnelt et al 
Gard ner et alj|200ll ; ICen et al. [ |2003l ; iNagamine et al 



2004 



20071 ). Only very recent works incorporat ed RT, primarily 



to account for the attenuation of the UVB dRazoumov et al 
20061: iPontzen et al.ll2008l;lFumagalli et al.ll201ll;lAltav et al 



2011 ; lMcQuinn et al.ll2011f ) and found a sharp transition be- 



tween optically thin and self-shielded gas that is expected 
from the exponential nature of extinction. 

The aforementioned studies focused mainly on redshifts 
z = 2 — 3, for which observational constraints are strongest, 
without investigating the evolution of the HI distribution. 
They found that the HI CDDF in current cosmological sim- 
ulations is in reasonable agreement with observations in a 
large range of HI column densities. Only at the highest HI 
column densities (i.e., Am > 10 21 cm -2 ) the agreement is 
poor. However, it is worth noting that the interpretation of 
these HI systems is complicated due to the complex physics 
of the ISM and ionization by local sources. Moreover, the ob- 
servational uncertainties are also larger for these rare high 
Am systems. 

In this paper, we investigate the cosmological HI distri- 
bution and its evolution during the last > 12 billion years 
(i.e., z < 5). For this purpose, we use a set of cosmologi- 
cal simulations which include star formation, feedback and 
metal-line cooling in the presence of the UVB. These simu- 
lations are based on th e Overwhelmingly L arge Simulations 
(OWLS) presented in ISchave et al.l |2010l 'l. To obtain the 
HI CDDF, we post-processed the simulations with RT, ac- 
counting for both ionizing UVB radiation and ionizing re- 
combination radiation (RR). In contrast to previous works, 
we account for the impact of recombination radiation explic- 
itly, by propagating RR photons. Using these simulations we 
study the evolution of the HI CDDF in the range of redshifts 
2 = — 5 for column densities Nui > 10 16 cm -2 . We discuss 
how the individual contributions from the UVB, RR and 
collisional ionization to the total ionization rate shape the 
HI CDDF and assess their relative importance at different 
redshifts. 

The structure of this paper is as follows. In [J2] we de- 
scribe the details of the hydrodynamical simulations and of 



the RT, including the treatment of the UVB and recombi- 
nation radiation. In |3]we present the simulated HI CDDFs 
and their evolution and compare them with observational 
constraints. In the same section we also discuss the contri- 
butions of different ionizing processes to the total ionization 
rate and provide fitting functions for the total photoioniza- 
tion rate as a function of density which reproduce the RT 
results. Finally, we conclude in [|4] 



2 SIMULATION TECHNIQUES 
2.1 Hydrodynamical simulations 

We use density fields from a set of cosmological simulations 
performed using a modified version of the smoothe d particle 
hydr odynamics code GADGET-3 (last described in lSpringell 
l2005h . The subgrid physics is identical to that used in the re f- 
erence simulation of the OWLS project l|Schave et alj|20"ich . 
Star formation is pressure depend ent and reproduces the 
obser ved Kennicutt-Schmidt law (|Schave fc Dalla Vecchial 
Chemical evolu tion is followed using the model of 



IWiersma et~aH (|2009bh . which traces the abundance evo- 
lution of eleven elements by following stellar evolution as- 
suming a IChabrierl (|2003T ) initial mass function. Moreover, 
a radiative heating an d cooling implementation based on 
W iersma et al.l (|2009aT ) calculates cooling rates element-by- 
element (i.e., using the above mentioned 11 elements) in the 
presence of the uniform c osmic microwave backgro und and 
the UVB model given by lHaardt fc Madaul l|200ll) . About 
40 per cent of the available kinetic energy in type II SNe 
is injected in winds with initial ve locity of 600 kms" 1 and 
a ma ss loading parameter rj = 2 l|Dalla Vecchia fc Schavel 
120081 ). 

We adopt fiducial cosmological parameters consis- 
tent with the most recent WMAP 7-year results: Q m = 
0.272, n h = 0.04 55, Qa = 0.728 a 8 = 0.81, n s = 0.967 
and h = 0.704 dKomatsu et all 120111) . We also use cos- 
mological simulations from the OWLS project which are 
performed with a cosmology consistent with WMAP 3-year 
values with O m = 0.238, Q, h = 0.0418, Q.a = 0.762, a s = 
0.74, n s — 0.951 and h — 0.73. We use those simulations to 
avoid expensive resimulation with a WMAP 7-year cosmol- 
ogy. Instead, we correct for the difference in the cosmological 
parameters as explained in Appendix iBl 

Our simulations have box sizes in the range L — 6.25 — 
100 comoving h~ Mpc and baryonic particle masses in the 
range 1.7 x 10 5 /i _1 M -8.7x 10 7 /i _1 M . The suite of sim- 
ulations allows us to study the dependence of our results on 
the box size and mass resolution. Characteristic parameters 
of the simulations are summarized in Table [1] 



2.2 Radiative transfer with TRAPHIC 

The RT is performed using TRAPHIC l|Pawlik fc Schavd 
I2OO8L l201ll ). TRAPHIC is an explicitly photon-conserving 
RT method designed to transport radiation directly on the 
irregular distribution of SPH particles using its full dynamic 
range. Moreover, by tracing photon packets inside a discrete 
number of cones, the computational cost of the RT becomes 
independent of the number of radiation sources. TRAPHIC is 
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Tabl e 1 ■ List of cosmological simulations used in this work. All the simulations use model ingredients identical to the reference simulation 
of ISchave et al From left to right the columns show: simulation identifier; comoving box size; number of dark matter particles 

(there are equally many baryonic particles); initial baryonic particle mass; dark matter particle mass; comoving (Plummer-equivalent) 
gravitational softening; maximum physical softening; final redshift; cosmology. The last column shows whether the simulation was post- 
processed with RT. In simulations without RT, the HI distribution is obtained by using a fit to the photoionization rates as a function 
of density measured from simulations with RT. 



Simulation 


L 


N 




"1dm 


e com 


e prop 


^end 


Cosmology 


RT 




(h- 1 Mpc) 




(fc^Ma) 




(/i^kpc) 




(h-'-kpc) 






L06N256 


6.25 


256 3 


1.7 x 10 5 


7.9 x 10 5 


0.98 


0.25 


2 


WMAP7 


/ 


L06N128 


6.25 


128 3 


1.4 x 10 6 


6.3 x 10 6 


1.95 


0.50 





WMAP7 


/ 


L12N256 


12.50 


256 3 


1.4 x 10 6 


6.3 x 10 6 


1.95 


0.50 


2 


WMAP7 


/ 


L25N512 


25.00 


512 3 


1.4 x 10 6 


6.3 x 10 6 


1.95 


0.50 


2 


WMAP7 




L06N128-W3 


6.25 


128 3 


1.4 x 10 6 


6.3 x 10 6 


1.95 


0.50 


2 


WMAP3 


/ 


L25N512-W3 


25.00 


512 3 


1.4 x 10 6 


6.3 x 10 6 


1.95 


0.50 


2 


WMAP3 




L25N128-W3 


25.00 


128 3 


8.7 x 10 7 


4.1 x 10 8 


7.81 


2.00 





WMAP3 


/ 


L50N256-W3 


50.00 


256 3 


8.7 x 10 7 


4.1 x 10 s 


7.81 


2.00 





WMAP3 


/ 


L50N512-W3 


50.00 


512 3 


1.1 x 10 7 


5.1 x 10 7 


3.91 


1.00 





WMAP3 


/ 


L100N512-W3 


100.00 


512 3 


8.7 x 10 7 


4.1 x 10 8 


7.81 


2.00 





WMAP3 





therefore particularly well-suited for RT calculation in cos- 
mological density fields with a large dynamical range in 
densities and large numbers of sources. In the following we 
briefly describe how TRAPHIC wor ks. More details, as well 
as various RT tests, can be found in iPawlik fc Schavd |2008, 

The photon transport in TRAPHIC proceeds in two 
steps: the isotropic emission of photon packets with a char- 
acteristic frequency v by source particles and their subse- 
quent directed propagation on the irregular distribution of 
SPH particles. The spatial resolution of the RT is set by the 
number of neighbors for which we generally use the same 
number of SPH neighbors used for the underlying hydrody- 
namical simulations, i.e., -/V ng b = 48. 

After source particles emit photon packets isotropi- 
cally to their neighbors, the photon packets travel along 
their propagation directions to other neighboring SPH parti- 
cles which are inside their transmission cones. Transmission 
cones are regular cones with opening solid angle Ati/Ntc 
and are centered on the propagation direction. The parame- 
ter A^tc sets the angular resolution of the RT, and we adopt 
Arc = 64. We demonstrate convergence of our results with 
the angular resolution in Appendix [C] Note that the trans- 
mission cones are defined locally at the transmitting particle, 
and hence the angular resolution of the RT is independent 
of the distance from the source. 

It can happen that transmission cones do not contain 
any neighboring SPH particles. In this case, additional parti- 
cles (virtual particles, ViPs) are placed inside the transmis- 
sion cones to accomplish the photon transport. The ViPs, 
which enable the particle-to-particle transport of photons 
along any direction independent of the spatially inhomoge- 
neous distribution of the particles, do not affect the SPH 
simulation and are deleted after the photon packets have 
been transferred. 

An important feature of the RT with TRAPHIC is the 
merging of photon packets which guarantees the indepen- 
dence of the computational cost from the number of sources. 
Different photon packets which are received by each SPH 
particle are binned based on their propagation directions in 
A^rc reception cones. Then, photon packets with identical 



frequencies that fall in the same reception cone are merged 
into a single photon packet with a new direction set by the 
weighted sum of the directions of the original photon pack- 
ets. Consequently, each SPH particle holds at most A^rc X N v 
photon packets, where N v is the number of frequency bins. 
We set A^rc = 8 for which our tests yield converged results. 

Photon packets are transported along their propaga- 
tion direction until they reach the distance they are al- 
lowed to travel within the RT time step by the finite 
speed of light, i.e., cAt. Photon packets that cross the 
simulation box boundaries are assumed to be lost from 
the computational domain. We use a time step At = 

1 ^ ( 6 . 25 \ b -iM P c ) (ife) (li|7) ; where A S ph is the 
number of SPH particles in each dimension. We verified that 
our results are insensitive to the exact value of the RT time 
step: values that are smaller or larger by a factor of two pro- 
duce essentially identical results. This is mostly because we 
evolve the ionization balance on smaller subcycling steps, 
and because we iterate for the equilibrium solution, as we 
discuss below. At the end of each time step the ionization 
states of the particles are updated based on the number of 
absorbed ionizing photons. 

The number of ionizing photons that are absorbed dur- 
ing the propagation of a photon packet from one particle to 
its neighbor is given by 5A4bs,i/ = &A/ln,i/[l — exp(— t(u))] 
where <5A/j n ,„ and r(y) are, respectively, the initial num- 
ber of ionizing photons in the photon packet with fre- 
quency v and the total optical depth of all the absorb- 
ing species. In this work we mainly consider hydrogen ion- 
ization, but in general the total optical depth is the sum 
T ( u ) ~ ~l2a Ta ( u ) °f the- optical depth of each absorbing 
species (i.e., a £ {HI, Hel, Hell}). Assuming that neighbor- 
ing SPH particles have similar densities, we approximate the 
optical depth of each species using r a {u) — cr a (v)n a dabs, 
where n a is the number density of species, d a b s is the ab- 
sorption distance between the SPH particle a nd its neigh- 
bor a nd <j a {v) is the absorption cross section l|Verner et al.l 
1996). Note that ViPs are deleted after each transmission, 
and hence the photons they absorb need to be distributed 
among their SPH neighbors. However, in order to decrease 
the amount of smoothing associated with this redistribution 
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of photons, ViPs are assigned only 5 (instead of 48) SPH 
neighbors. We demonstrate convergence of our results with 
the number of ViP neighbors in Appendix [C] 

At the end of each RT time step, every SPH particle has 
a total number of ionizing photons that have been absorbed 
by each species, AA/kbs.^^). This number is used in order to 
calculate the photoionization rate of every species for that 
SPH particle. For instance, the hydrogen photoionization 
rate is given by: 

Fm = ^HiA/nAt ' (1) 

where Ah is the total number of hydrogen atoms inside the 
SPH particle and ?7hi = n H i/ n H i s the hydrogen neutral 
fraction. 

Once the photoionization rate is known, the evolution of 
the ionization states is calculated. For instance, the equation 
which governs the ionization state of hydrogen is 



— =a ml n B (l- 



tjhi) — ^Hi(r H i + r e ,Hn e ), 



(2) 
is the 



where n e is the free electron number density, r e: H 
collisional ionization rate and ami is the HII recombination 
rate. The differential equations which govern the ionization 
balance (e.g., equation^) are solved using a subcycling time 

/ (T ion + Tree J , 



' ion ' rcc 



step, St = min(/T eq ,At) where r e , 
and / is a dimensionless factor which controls the integration 
accuracy (we set it to 10 -3 ), r rec = 1/ ^2 i n e at and Tio» = 
1 / yV (Tj + n e r e ,i). The subcycling scheme allows the RT 
time step to be chosen independently of the photoionization 
and recombination time scales without compromising the 
accuracy of the ionization state calculation^]. 

We employ separate frequency bins to transport UVB 
and RR photons. Because the propagation directions of pho- 
tons in different frequency bins are merged separately, this 
allows us to track the individual radiation components, i.e., 
UVB and RR, and to compute their contributions to the 
total photoionization rate. The implementation of the UVB 
and the recombination radiation is described in § 12.31 and 
§H3] below. 

At the start of the RT, the hydrogen is assumed 
to be ne utral. In addition, we use a common simplifica- 



tion ( e.g., Fauchcr-Gigucrc et alj |2009: Mc Quinn fc Switzerl 



l20ld : lAltav et alj|201lf i ibv assuming a hydrogen mass frac- 
tion of unity, i.e., we ignore helium (only for the RT). To 
calculate recombination and collisional ionizations rates, we 
set, in post-processing, the temperatures of star- forming gas 
particles with densities n H > 0.1cm" ~ 3 to Tism = 10 4 K, 
which is typical of the observed warm-neutral phase of the 
ISM. This is needed because in our hydrodynamical sim- 
ulations the star-forming gas particles follow a polytropic 
equation of state which defines their effective temperatures. 
These temperatures are only a measure of the imposed 
pressure and do not represent physical temperatures (see 



2 Other considerations prevent the use of arbitrarily large RT 
time steps. The RT assumes that species fractions and hence 
opacities do not evolve within a RT time step. This approximation 
becomes increasingly inaccurate with increasing RT time steps. 
Note that in this work, we iterate for ionization equilibrium which 
help to render our results robust against changes in the RT time 
step, as our convergence studies confirm. 



ISchave fc Dalla Vecchiafe oOS ) . To speed up convergence, the 
hydrogen at low densities (i.e., n H < 10 -3 cm -3 ) or high 
temperatures (i.e., T > 10 5 K) is assumed to be in ioniza- 
tion equilibrium with the UVB and the collisional ionization 
rate (see Appendix IA2|) . Typically, the neutral fraction of 
the box and the resulting HI CDDF do not evolve after 2-3 
light-crossing times (the light-crossing time for the extended 
box with Lbox = 6.25 comoving fa -1 Mpc is w 7.5 Myr at z 
= 3). 



2.3 Ionizing background radiation 

We simulate the ionizing background radiation as plane- 
parallel radiation entering the simulation box from its sides. 
At the beginning of each RT step, we generate a large num- 
ber of photon packets, Nbg, on the nodes of a regular grid 
at each side of the simulation box and set their propaga- 
tion directions perpendicular to the sides. The number of 
photon packets is chosen to obtain converged results. Fur- 
thermore, to avoid numerical artifacts close to the edges of 
the box, we use the periodicity of our simulations to extend 
the simulation box by the typical size of the region where we 
generate the background radiation (i.e., 2% of the box size 
from each side). These extended regions are excluded from 
the analysis, thereby removing the artifacts without losing 
any information contained in the original simulation box. 

The photon content of each packet is normalized such 
that in the absence of any absorption (i.e., assuming the 
optically thin limit), the total photon density of the box 
corresponds to the desired uniform hydrogen photoioniza- 
tion rate. If we assume that all the photons with frequencies 
higher than i^Heii are absorbed by helium, then the hydrogen 
photoionization rate can be written as: 



"Hell j 

4n — rJHi,i/ dv 



4n a 



"HI 



h 



"Hell J 

— dv, 
v 



(3) 



"HI 



where J„ is the radiation intensity (in units 
erg cm - s~ sr~ Hz - ), vm and i/emi are respectively 
the frequency at the Lyman-limit and the frequency at the 
Hell ionization edge, and cthi.i/ is the neutral hydrogen 
absorption cross-section for ionizing photons. In the last 
equation we have defined the gray absorption cross-section, 



j "hi ; i 

r H ° n Ju/v dv • 

J "HI ' 



(4) 



The radiation intensity is related to the photon energy 
density, w„, 

T u v c n„ hv c 
J " = —a — = — ; . ( 5 ) 

47T 47T 

where n u is the number density of photons inside the box. 
Combining Equations 13151 yields 



Tuvb 



(6) 



where n„ HI is the number density of ionizing photons inside 
the box. The total number of ionizing photons in the box is 
therefore given by 

-£-box 



£hox = ru, 6 TVb 



"HI -^box 



At' 



(7) 
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Table 2. Hydrogen photoionization rate, absorption cross-s ection, equivalent gray ap proximation frequency and the self-shielding density 
threshold (i. e.. based on equationlT3ll for th ree UVB models: lHaardt fc Madaul l|200lf ) (HM01; used in this work) . lHaardt fc Madaul ((20121 ) 
(HM12) and lFaucher-Giguere et al.l {2009) (FG09) at different redshifts. For the calculation of the photoionization rate and absorption 
cross-sections, only photons with energies below 54.4 eV are taken into account, effectively assuming that more energetic photons are 
absorbed by He. 



Redshift 


UVB 


Tuvb (s 


X ) 


°Vhi ( cm2 ) 


S eq (eV) 


™H,sSh (cm 3 ) 


z = 


HM01 


8.34 x 10" 


14 


3.27 x 10- 18 


16.9 


1.1 x lO" 3 




tTM1 9 


9 97 v 1 O - 


14 


Z.Oo X 1U 


lO 1 

lo. 1 


O.l X 1U 




FG09 


3.99 x 10" 


14 


2.59 x 10~ 18 


18.3 


7.7 x 10" 4 


2 = 1 


HM01 


7.39 x 10" 


13 


2.76 x 10" 18 


17.9 


5.1 x 10" 3 




HM12 


Q AO v 1(1" 


13 


9 fi9 v i n — is 


18 2 


O.O X _1_U 




FG09 


3.03 x 10" 


13 


2.37 x 10" 18 


18.8 


3.1 x 10" 3 


2 = 2 


HM01 


1.50 x 10" 


12 


2.55 x 10 -18 


18.3 


8.7 x 10" 3 




HMf 2 


fi Q2 y 1(1- 

o.yo x 


13 


Q fil y i n — 18 


18.2 


6.1 x 10 — 3 




FG09 


6.00 x 10- 


13 


2.27 x 10" 18 


19.1 


5.1 x lO" 3 


2 = 3 


HM01 


1.16 x 10" 


12 


2.49 x 10" 18 


18.5 


7.4 x 10" 3 




HM12 


8.74 x 10- 


13 


2.61 x 10" 18 


18.2 


6.0 x 10~ 3 




FG09 


5.53 x 10- 


13 


2.15 x 10- 18 


19.5 


5.0 x lO" 3 


2 = 4 


HM01 


7.92 x 10- 


13 


2.45 x 10" 18 


18.6 


5.8 x 10" 3 




HM12 


6.14 x 10- 


13 


2.60 x 10" 18 


18.3 


4.7 x 10- 3 




FG09 


4.31 x 10- 


13 


2.02 x 10" 18 


19.9 


4.4 x lO" 3 


2 = 5 


HM01 


5.43 x 10- 


13 


2.45 x 10~ 18 


18.6 


4.5 x 10- 3 




HM12 


4.57 x 10- 


13 


2.58 x 10" 18 


18.3 


3.9 x lO" 3 




FG09 


3.52 x 10- 


13 


1.94 x 10" 18 


20.1 


4.0 x 10" 3 



where n 7 is the number of ionizing photons carried by each 
photon packet. Now we can calculate the photon content of 
each packet that must be injected into the box during each 
step in order to achieve the desired HI photoionization rate: 



1 UVB 



6 ov HI N\ 



(8) 



We use the redsh ift-dependent UVB spectrum of 

The 



Haardt & Madau (2001) to calculate Tuvb and a u 



Haardt fc Madaul ((2001) UVB model successfully repro- 
duces the relative strengths of the o bserved metal absor ption 
lines in the intergalactic medium (Aguirrc ct al. 2008j) and 
has been used to calculate heating/cooling in our cosmolog- 
ical simulations 0. 

To reduce the computational cost, we treat the multi- 
frequency problem in the gray approximation. In other 
words, we transport the UVB radiation using a single fre- 
quency bin, inside which photons are absorbed using the 
gray cross-section a VHI defined in equation [3] Note that the 
gray approximation ignores the spectral hardening of the ra- 
diation field that would occur in multifrequency simulations. 
In Appendix [D] we confirm the validity of our assumptions 
by repeating one of our simulations using multiple frequency 



3 Note that during the hydrodynamical simulations, photoheat- 
ing from the UVB is applied to all gas particles. This ignores 
the self-shielding of hydrogen atoms against the UVB that occurs 



at densities nn > 10 



10 cm 3 . This inconsistency, which 



could affect both collisional ionization rates and the small-scale 
structure of the absorbers, has been foun d to have no signifi- 
cant impact on the simulated HI CDDF l lPontzen et al.l 120081 ; 
iMcQuinn fc SwitzedliuTf] ; lAltav et al.1l201ll) . 



bins, and also explicitly accounting for the absorption of 
photons by helium. 

Hydrogen photoionization rates and average absorp- 
tion cross-sections for UVB radiation at different redshifts 
ar e listed in Table [Z] for our fiducial UVB model based 
onlHaardt fc Madauf(|200ir i together with lHaardt fc Madaul 
i|20l2l) . The photoionization rate peaks at 2 ~ 2 — 3 in 
both models and the equivalent effective photon energ^Q of 
the background radiation changes only weakly with redshift, 
compared to the total photoionization rate. 



2.4 Recombination radiation 

Photons produced by the recombination of positive ions and 
electrons can also ionize the gas. If the recombining gas is 
optically thin, recombination radiation can escape and its 
ionizing effects can be ignored (i.e., the so-called Case A). 
However, for regions in which the gas is optically thick, the 
proper approximation is to assume the ionizing recombina- 
tion radiation is absorbed on the spot. In this case, the 
effective recombination rate can be approximated by ex- 
cluding the transitions that produce ionizing photons (e.g., 
lOsterbrock fc Ferlandll2006l ). This scenario is usually called 
Case B. A possible way to take into account the effect of 
recombination radiation is to use Case A recombination at 
low densities and Case B recombination at high densities 
fe.g.. lAltav et al.ll201ll ; IMcQuinn et al.ll201ll ). but this will 
be inaccurate in the transition regime. 

4 We defined the equivalent effective photon energy, £ cq , which 
corresponds to the absorption cross section, av HI , as: £ eq = 



13.6 eV 



-1/3 



where ao = 6.3 X 10 18 cm 2 . 
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In this work we explicitly treat the ionizing photons 
emitted by recombining hydrogen atoms and follow their 
propagation through the simulation box. This is facili- 
tated by the fact that the computational cost of RT with 
TRAPHIC is independent of the number of sources. This 
is particularly important noting that every SPH particle is 
potentially a source. The photon production rates of SPH 
particles depend on their recombination rates and the radi- 
ation is emitted isotropically once at the beginning of every 
RT time step (see Raicevic et al. in prep, for full details). 

We do not take into account the redshifting of the re- 
combination photons by peculiar velocities of the emitters, 
or the Hubble flow. Instead, we assume that all recombina- 
tion photons are monochromatic with energy 13.6 eV. 

2.5 The HI column density distribution function 

In order to compare the simulation results with observa- 
tions, we compute the CDDF of neutral hydrogen, f(Nm, z), 
which is defined as the number of absorbers per unit column 
density, per unit absorption length, dX: 



/(AW) 



d 2 n 



d 2 n H{z) 



1 



dNmdX dNmdz Ho (1 + zf 



(9) 



We project the HI content of the simulation box along 
each axis onto a grid with 5000 2 or 10000 2 pixels (for the 
128 3 - 256 3 and 512 3 simulations, respectively )Q. This is 
done using the actual kernels of SPH particles and for each 
of the three axes. The projection may merge distinct systems 
along the line of sight. However, for the small box sizes and 
high column densities with Ahi > 10 17 cm -2 , which are the 
focus of this work, the chance of overlap between multiple 
absorbing systems in projection is negligible. Based on our 
numerical experiments, we expect that this projection ef- 
fect starts to appear only at Ahi < 10 16 cm -2 if one uses 
a single slice for the projection of the entire L50N512-W3 
simulation box at z — 3. To make sure our results are insen- 
sitive to this effect, we use, depending on redshift, 25 or 50 
slices for projecting the L50N512-W3 simulation. 

To produce a converged /(iVm, z) from simulations, one 
needs to use cosmological boxes that are large enough to 
capture the relevant range of over-densities. This is partic- 
ularly demanding at very high H I column densities: for in- 
stance, Ivan de Voort et al.l (|2012l ') showed that most of the 
gas with Am > 10 21 cm~ 2 resides in galaxies with halo 
masses > 10 11 Mq which are relatively rare. As we show in 
Appendix(B] the box size required to produce a converged HI 
CDDF up to Am ~ 10 22 cm" 2 is L > 50 comoving /i _1 Mpc. 
Simulating RT in such a large volume is expensive. However, 
as we show in H3.41 at a given redshift the photoionization 
rates are fit very well by a function of the hydrogen num- 
ber density. This relation is conserved with respect to both 
box siz^fl and resolution and can therefore be applied to our 
highest resolution simulation (i.e., L50N512-W3), allowing 

5 Using 5000 2 cells, the corresponding cell size is similar to the 
minimum smoothing length, and ~ 100 times smaller than the 
mean smoothing length, of SPH particles at z = 3 in the L06N128 
simulation. 

6 One should note that the box size can indirectly change the 
resulting photoionization rate profile. For instance, self-shielding 
can be affected by collisional ionization, which become stronger at 



us to keep the numerical cost tractable. Finally, since re- 
peating the high-resolution simulations is expensive, we ap- 
ply a redshift-independent correction which accounts for the 
difference between the WMAP year 3 parameters used for 
L50N512-W3 simulation and the WMAP year 7 values. This 
is done by multiplying all the HI CDDFs produced based on 
the WMAP year 3 cosmology by the ratio between the HI 
CDDFs for L25N512 and L25N512-W3 at z = 3. 



2.6 Dust and molecular hydrogen 

Dust extinction is a physical processes not treated in our 
simulations. Assuming a constant dust-to-gas ratio, the typ- 
ical dust absorption cross-section per atom is orders of mag- 
nitudes lower than the typic al hydrogen absorption cross - 
section for ionizing photons ( Wei ngartner fc Drainel [2001 1. 
In other words, the absorption of ionizing photons by dust 
particles is not significant compared to the absorption by 
the neutral hydrogen. Consequently, as also found in cos- 
molo gical simulations with ionizing radiation (|Gnedin et al.l 
2008), dust absorption does not noticeably alter the overall 
distribution of ionizing photons and hydrogen neutral frac- 
tions. 

The observed cut-off in the abundance of very high 
A'hi systems may be rela ted to the con v ersion of atomic 
hydrogen into H2 (e.g., Schavel 2001bl; Krumholz et al.l 



l2009bl;|Prochas-ka fc Wolfd 20091: lAltav et alj|201ll ). Follow- 
ing I Altav et al.l (|201lf ) and lDuffv et al.l (|2012l '). we adopt an 
observationally driven scaling rel ation between gas pressure 
and hydrogen molecular fraction l|Blitz fc Rosolowsky||2006h 
in post-processing, which reduces the amount of observable 
HI at high densities. This scaling relation is based on obser- 
vations of low-redshift galaxies and may not cover the low 
metallicities relevant for higher redshifts. This could be an 
issue, since the HI-H2 relation is known to be sens itive to 
the dust content and hence to the m etallicity (e.g.. ISchavd 
l2001bl . 12004 iKrumholz et~alll2009ah . 



3 RESULTS 

In this section we report our findings based on various RT 
simulations which include UVB ionizing radiation and dif- 
fuse recombination radiation from ionized gas. As we demon- 
strate in H3.3I the dependence of the photoionization rate on 
density obtained from our RT simulations shows a generic 
trend for different resolutions and box sizes. Therefore, we 
can use the results of RT calculations obtained from smaller 
boxes (e.g., L06N128 or L06N256) which are computation- 
ally cheaper, to calculate the neutral hydrogen distribution 
in larger boxes. The last column of Table [T] indicates for 
which simulations this was done. 

In the following, we will first present the predicted HI 
CDDF and compare it with observations. Next we discuss 
other aspects of our RT results and the effects of ionization 
by the UVB, recombination radiation and collisional ioniza- 
tion on the resulting HI distributions at different redshifts. 



lower redshifts and whose importance depends on the abundance 
of massive objects, which is more sensitive to the box size. 
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Figure 1. CDDF of neutral gas at different redshifts in the presence of the UVB and diffuse recombination radiation for L50N512-W3. 
A column density dependent amplitude correction has been applied to make the results consistent with WMAP year 7 cosmological 
parameters. The observati onal data points rep resent a compilation of various quas a r abso rption line observati ons at high redshif t s (i.e. , 
z = [1.7,5.5]) taken f rom IPeroux et alj J2005h w ith z = [1.8.3.5]. lo'Meara et all l|2007T l with z = [1.7. 4.5], iNoterdaeme et al.1 l|2009h 
with z = [2 2,5.5] and iProchaska fc Wolfe! J2009I ) with z = [2.2,5.5]. The colored data points in the top-left corner of the left panel are 
taken from lKim et al with z = [2.9,3.5] and z = [1.7,2.4] for the yello w crosses and oran ge diamonds, respectively. The orange 

filled circles show the best-fit based on the low-redshift 21-cm observations of IZwaan et al.l ll2005h . The high column density end of the 
HI distribution is magnified in the right panel and for clarity only the simulated HI CDDF of redshifts z = 1, 3 & 5 are shown. The 
top-section of each panel shows the ratio between the HI CDDFs at different redshifts and the HI CDDF at z = 3. The simulation results 
are in reasonably good agreement with the observations and, like the observations, show only a weak evolution for Lyman Limit and 
weak damped Lyct systems below 2 = 3. 



3.1 Comparison with observations 

In Figure [1] we compare the simulation results with a com- 
pilation of observed HI CDDFs, after converting both to the 
WMAP year 7 cosmology. The data points with error bars 
show results from high-redshift (z = 1.7 — 5.5) QSO absorp- 
tion line studies and the orange fil l ed cir cles show the fitting 
function reported bv lZwaan et al] l|2005h based on 21-cm ob- 
servations of nearby galaxies. The latter observations only 
probe column densities Am > 10 19 cm" 2 . 

We note that the OWLS simula tions have a l ready been 
shown to agree with observations bv lAltav et"aH (|201ll ). but 
only for z — 3 and based on a different RT method (see Ap- 
pendix [C3l for a comparison). Overall, our RT results are also 
in good agreement with the observations. At high column 
densities (i.e., Am > 10 17 cm -2 ) the observations probing 
< z < 5.5 are consistent with each other. This implies 
weak or no evolution with redshift. The simulation is con- 
sistent with this remarkable observational result, predicting 
only weak evolution for 10 17 cm~ 2 < Am < 10 21 cm -2 (i.e., 



Lyman limit systems, LLSs, and weak Damped Ly-o sys- 
tems, DLAs) especially at z < 3. 

The simulation predicts some variation with redshift for 
strong DLAs (Am > 10 21 cm -2 ). The abundance of strong 
DLAs in the simulations follows a similar redshift- dependent 
trend as the average star formati on density in our s imula- 
tions which peaks at z « 2 — 3 l|Schave et al]|201ol ). This 
result is consistent with the DLA evolution found by ICenl 
(2012) in two zoomed simulations of a cluster and a void. 
One should, however, note that at very high column densi- 
ties (e.g., A^hi > 10 21 ' 5 cm -2 ) both observations and simu- 
lations are limited by small number statistics and the sim- 
ulation results are more sensitive to the adopted feedback 
scheme (Altay et al. in prep.). Moreover, as we will show 
in Rahmati et al. (in prep.), including local stellar ionizing 
radiation can decrease the HI CDDF by up to « 1 dex for 
Nm > 10 21 cm -2 , especially at redshifts z ~ 2 — 3 for which 
the average star formation activity of the Universe is near 
its peak. 
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At low column densities (i.e., iVm < 10 17 cm -2 ) the 
simulation results agree very well with the observations. 
This is apparent from the agreement between the simulated 
f(Nui,z) at z = 3 and 2 = 4, and the o bserved values for 
redshifts 2.9 < z < 3.5 jKim et al.l 120021 ) which are shown 
by the yellow crosses in the left panel of Figure [T] The sim- 
ulated f(Nm,z) at lower and higher redshifts deviate from 
those at z as 3 showing the abundance of those systems de- 
creases with decreasing redshift and remains nearly constant 
at z < 2. This is co nsistent with the Ly-q forest observations 
at lower redshifts dKim et alj 120021; jjanknecht et all l200rj; 
Lehner et ai1l2007l : jpTochaska fc Wolfdl2009l : iRibaudo et al l 
20111 ). as illustrated with the orange diamonds which cor- 
respond to z w 2 observations, in the top-left corner of the 
left panel in Figure [1] 

The evolution of the HI CDDF with redshift results 
from a combination of the expanding Universe and the 
growing intensity of the UVB radiation down to redshifts 
z « 2 — 3. At low redshifts (i.e., z « 0) the intensity of 
the UVB radiation has dropped by more than one order of 
magnitude leading to higher hydrogen neutral fractions and 
higher HI column densities. However, as we show in £|3 . 51 at 
lower redshifts an increasing fraction of low-density gas is 
shock-heated to temperatures sufficiently high to become 
collisionally ionized and this compensates for the weaker 
UVB radiation at low redshifts. 

The simulated HI CDDFs at all redshifts are consis- 
tent with each other and the observations. However, as il- 
lustrated in the right panel of Figure [T] there is a as 0.2 dex 
difference between the simulation results and the observa- 
tions of LLS and DLAs at all redshifts. We found that the 
normalization of the HI CDDF in those regimes is sensitive 
to the adopted cosmological parameters (see Appendix |B)| . 
Notably, the cosmology consistent with the WMAP 7 year 
results that is shown here, produces a better match to the 
observations than a cosmology based on the WMAP 3 year 
results with smaller values for fib and ag . This suggests that 
a higher value of as may explain the small discrepancy be- 
tween the simulation results and the observations. 



3.2 The shape of the HI CDDF 

The shape of the HI CDDF is determined by the distri- 
bution of hydrogen and by the different ionizing processes 
that set the hydrogen neutral fractions of the absorbers. 
One can assume that over-dense hydrogen resides in self- 
gravitating systems that are in local hydrostatic equilibrium. 
Then, the typical scales of the systems can be calculated as 
a functi on of the gas density based on a Jeans scaling ar- 
gument l|Schavell2001al ). Assuming that absorbers have uni- 
versal baryon fractions (i.e., f s = j^-) and typical temper- 
atures of Ti = (T/10 4 K) ~ 1 (i.e., collisional ionization is 
unimpor tant), one can calculate the total hydrogen column 
density (|Schavdl2001ah : 



i|Schavell200lah : 



N H ~ 1.6 x 10 21 cm -2 nj 2 



T 



1/2 



0.17 



1/2 



(10) 



Assuming that the gas is highly ionized and in ionization 
equilibrium with the ambient ionizing radiation field with 
the photoionization rate, F_i2 = T/10 -12 s _1 , one gets 



Nm ~ 2.3 x 10 13 cm -2 



T -0.26 r -i 
X J 4 1 _ 12 



10- 5 cm" 3 

1/2 



3/2 



0.17 



(11) 



At high densities where the gas is nearly neutral, equa- 
tion (|10[) provides a relation between JVhi and n H . Equation 
(|11[) on the other hand, gives the relation for optically thin, 
highly ionized gas. The latter is derived assuming that the 
UVB photoionization is the dominant source of ionization, 
which is a good assumption at high redshifts and explains 
the relation between d ensity and column density in Lyq 



forest simulations (e. g., Dave et "all 20ld; Attav et ai]|201ll ; 



iMcQuinn et all l201ll ; iTepper-Garcia et alj|2012h . However, 
as we will show in the following sections, photoionization 
domination breaks down at lower redshifts where collisional 
ionization plays a significant role. 

The column density at which hydrogen starts to be 
self-shielded against the UVB radiation follows from setting 
Thi = 1: 



N: 



HI SSh 



■4 x 10 17 cm" 2 



2.49 x 10- 18 cm- 



(12) 



which can be used together with equation (|11[) to find the 
typical densities at w hich the self-shielding begins (e.g., 
iFurlanetto et alj|2005ft : 



n H ssh ~ 6.73 x 10~" 3 cm~ 3 



x T. 



,0.17 t-,2/3 



0.17 



ff »m 

2.49 x 10" 18 cm 2 

1/3 



-2/3 



(13) 



These relations are compared with the n HI -weighted to- 
tal hydrogen number density as a function of JVhi in the 
L06N256 simulation at z — 3 in the left panel of Figure [2] 
The solid curve shows the median and the red (blue) shaded 
area represents the central 70% (90%) percentile. The di- 
agonal gray solid line which converges with the simulation 
results at low column densities, shows equation (|11[) and the 
steeper gray dotted line which converges with the simulation 
results at high column densities is based on equation (|10[) . 
The agreement between the expected slopes of the n H — Nm 
relation and the simulations at low and high column densi- 
ties confirms our initial assumption that hydrogen resides in 
self-gravitating systems which are close to local hydrostatic 
equilibriunfl 

As expected from equation (|13[) . at low densities the gas 
is optically thin and follows the Jeans scaling relation of the 
highly ionized gas. At n H > 0.01 cm -3 however, the relation 
between density and column density starts to deviate from 
equation (|11[) and approaches that of a nearly neutral gas. 



7 One should note that the above mentioned Jeans argument 
provides an order of magnitude calculation due to its simplify- 
ing assumptions (e.g., uniform density, universal baryon fraction, 
etc.). Although we may expect the predicted scaling relations to 
be correct, the very close agreement of the normalization with the 
simulations at low densities is coincidental. As the steeper gray 
dotted line which is based on equation (1 1 P shows, the simulated 
TVhi for a given ri H is pb 0.5 dex higher than implied by the Jeans 
scaling for the nearly neutral case (i.e., steep, gray solid line). 
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Figure 2. Left: n HI -weighted total hydrogen number density as a function of Nm- The brown solid curve shows the RT results and 
the purple dotted curve shows the optically thin limit. Blue dot-dashed, dot-dot-dot-dashed and long dashed curves assume models with 
self-shielding density thresholds of nn,SSh = 10" 1 , 10" 2 and 10 — 3 cm - 3 , respectively. All of the above mentioned curves show the 
median of the n HI -weighted total hydrogen number density at a given Nm . The gray thin lines show the expected Jeans scaling relations 
for optically thin gas (equation 1101 diagonal solid line) and for neutral gas (equation 1111 steeper dotted line). A second solid line with 
the same slope expected from equation JTTJ but a different normalization is illustrated by the second solid line which is identical to 
the dotted line but shifted by 0.5 dex to higher TVhi- The pink and blue shaded areas in the right panel indicate the 70% and 99% 
scatter respectively, while the solid curves shows the median for the RT result. All the other curves are also medians. This shows that 
the n H — Nm relationship can be explained by the Jeans scaling and that the flattening in the CDDF is due to self-shielding. Bight: HI 
CDDF in the presence of the UVB and diffuse recombination radiation for simulation L06N256. Simulations shown with different curves 
are identical to the left panel. In addition, the effect of H2 formation is shown by the green solid curve which deviates from the brown 
solid curve at Nm > 3 X 10 21 cm -2 . Finally the red dashed curve, which is indistinguishable from the brown solid curve, shows the result 
of assuming the median of the photoionization rate profile of the RT results to calculate the neutral fractions (see £13.31 and Appendix 
lAH , The top-section in the right panel shows the ratio between different HI CDDFs and the one resulting from the RT simulations. 



Consequently, for densities above the self-shielding thresh- 
old the HI column density increases rapidly over a narrow 
range of densities, leading to a flattening in the n H — Nui 
relation and in the resulting HI CDDF at Nm > 10 18 cm" 2 
(see Figure [2}. The results from the RT simulation deviate 
from the magenta dotted lines, which are obtained assuming 
optically thin gas, at Nm > 4 x 10 17 cm~ 2 . As the dotted 
line in the right panel of Figure [2] shows, in the absence 
of self-shielding, the slope of f(Nm,z) oc Nm* 3 is constant 
all the way up to DLAs at fHhya ~ —1.6. However, because 
of self-shielding, the HI CDDF flattens to /3lls ~ —1.1 at 
10 18 cm" 2 < iV H i < 10 20 cm" 2 in the RT simulation (solid 
curve). These predicted slopes are in excellent agreement 
with the latest observational constraints of /?Lya > —1.6 
for 10 15 cm" 2 < A^hi < 10 17 cm" 2 to /3 LLS « -1 in 
the LLS regime (|Q'Meara et al.ll2012j ). We also note that 
Phya > —1.6 is predicted to be almost the same for all red- 
shifts . which agrees well with observations Jjanknecht et all 
120061 ; iLehner et al1l2007l ; iRibaudo et al.ll201ll ). 

At densities n H > 0.1 cm -3 the gas is nearly neu- 



tral and the Jeans scaling in equation (|1UI) controls the 
n H — Nm relation. Consequently, the rate at which Nm 
responds to changes in n H slows down, causing a steep- 
ening in the resulting f(Nm,z) in the DLA range (i.e., 



However, as the thick solid curve 



in the right panel of Figure [2] illustrates, the slope of 
f(Nm,z) remains constant for Am = 10 21 — 10 22 cm -2 . 
This is in contrast with observed tr ends indicating a sharp 



cut-off at Am > 3 x 10 21 c m 2 JProchaska et al 



2010; 



lO'Meara et all |2012| ; but see iNoterdaeme et al.l I2012T ). At 



those column densities a large fraction of hydrogen is ex- 
pected t o form H2 molecules and be absent from HI obser- 
vations |Schavdl2001bl ; iKrumhob et al. l l2009d ; [Aitav et al.l 
l201ll ). As the thin solid line in the right panel of Figure 
[2] shows, accounting for H2 using the empirical relation be- 
tween H2 fraction and pres sure, based on z = observations 
|Blitz fc Rosolowsky||2006l ). does reprod uce a sharp cut-off. 
If the observed relation does not cut off (|Noterdaeme et alj 
2012), then this may imply that H2 fractions are lower at 
z — 3 than at z = 0. We also note that the ionizing effect 
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Figure 3. The hydrogen neutral fraction {left) and the photoionization rate (right) as a function of hydrogen number density do not 
change by varying the simulation box size or mass resolution. This is shown for different simulations at z = 3 in the presence of the UVB 
and recombination radiation. Purple solid, blue dashed and red dot-dashed lines show, respectively, the results for L12N256, L06N128 
and L06N256. The green dotted line indicates the results for the L06N128 simulation if the gas is assumed to be optically thin to 
the UVB radiation (i.e., no RT calculation is performed). The deviation between the optically thin hydrogen neutral fractions and RT 
results at n H > 10 — 2 cm -3 shows the impact of self-shielding. The lines show the medians and the shaded areas indicate the 15% — 85% 
percentiles. At the top of each panel we show HI column densities corresponding to each density. 



of local sources (Rahmati et al. in prep.), increasing the ef- 
ficiency of stellar feedback, e.g., by using a top-heavy IMF, 
and AGN feedback can also affect these high HI column 
densities (Altay et al. in prep.). 

To first order, one can mimic the effect of RT by as- 
suming gas with n H < rm.ssh to be optically thin (i.e., 
Case A recombination) and gas with n H > nn,ssh to be 
fully neutral. Simulations with three different self-shielding 
density thresholds are shown in Figure The dot-dashed, 
dot-dot-dot-dashed and long dashed curves correspond to 
«H,ssh = 10 _1 , 10 -2 and 10 -3 cm -3 , respectively. Although 
all of these simulations predict the flattening of f(Nm,z), 
they produce a transition between optically thin and neu- 
tral gas that is too steep. In contrast, the RT results show 
a transition between highly ionized and highly neutral gas 
that is more gradual, as observed. 



3.3 Photoionization rate as a function of density 

Figure [3] illustrates the RT results for neutral fractions and 
photoionization rates as a function of density in the pres- 
ence of UVB radiation and diffuse recombination radiation 
for the L06N128, L06N256 and L12N256 simulations at z 
= 3. For comparison, the results for the optically thin limit 
are shown by the green dotted curves. The sharp transition 
between highly ionized and neutral gas and its deviation 
from the optically thin case are evident in the left panel. 
This transition can also be seen in the photoionization rate 
(right panel) which drops at n H > 0.0 1 cm" 3 , consistent 
with e quation (|13l> and previous studies dTaiiri fc Umemural 
1998; Razoum ov et al.l 20061; iFau chcr- Gigu ere et al.l 2010; 
Nagamine et al.l |2010| ; iFumagalli et all 1201 it lAltav et all 
201ll ). 

The medians and the scatter around them are insensi- 
tive to the resolution of the underlying simulation and to the 
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Figure 4. Ionization rates due to different sources of ioniza- 
tion as a function of hydrogen number density. Blue solid, green 
dashed and red dotted curves show, respectively, the UVB pho- 
toionization rate, the recombination radiation photoionization 
rate and the collisional ionization rate. The curves show the me- 
dians and the shaded areas around the medians indicate the 
15% — 85% percentiles. HI column densities corresponding to 
each density are shown along the top x-axis. While the UVB is 
the dominant source of ionization below the self-shielding (i.e., 
n H < 10 -2 cm — 3 ), recombination radiation dominates the ion- 
ization at higher densities. 



box size. This suggests that one can use the photoionization 
rate profile obtained from the RT simulations for calculat- 
ing the hydrogen neutral fractions in other simulations for 
which no RT has been performed. 

Moreover, as we show in £|3 . 5 1 the total photoionization 
rate as a function of the hydrogen number density has the 
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same shape at different redshifts. This shape can be charac- 
terized by three features: i) a knee at densities around the 
self-shielding density threshold, ii) a relatively steep fall-off 
at densities higher than the self-shielding threshold and iii) 
a flattening in the fall-off after the photoionization rate has 
dropped by ~ 2 dex from its maximum value which is caused 
by the RR photoionization. These features are captured by 
the following fitting formula: 



= 0.98 



1 + 



( n H 



V n H,S 



+0.02 



1 + - 



«-H,SSh 



(14) 



where Fuvb is the background photoionization rate and 
Tphot is the total photoionization rate. Moreover, the self- 
shielding density threshold, nH,ssh, is given by equation (|13[) 
and is thus a function of Tuvb and o"„ HI which vary with 
redshift. As explained in more detail in Appendix IA1I the 
numerical parameters representing the shape of the pro- 
file are chosen to provide a redshift independent best fit 
to our RT results. In addition, the parametrization is based 
on the main RT related quantities, namely the intensity of 
UVB radiation and its spectral s hape. It can therefore be 
used for UVB models similar to the Haardt & Madau (2001) 
model we used in this work (e.g., iFaucher-Giguere et al.l 
120091 : lHaardt fc Madaull2012l 'l. For a given UVB model, one 
only needs to know Fuvb and av HI in order to determine the 
corresponding nn,ssh from (| 13p (see also Table[5]). Then, af- 
ter using equation (|14|l to calculate the photoionization rate 
as a function of density, the equilibrium hydrogen neutral 
fraction for different densities, temperatures and redshifts 
can be readily calculated as explained in Appendix I All 

We note that our fitting function is only accurate for 
photoionization dominated cases. As we show in H3.51 at 
z ~ the collisional ionization rate is greater than the 
total photoionization rate around the self-shielding density 
threshold. Consequently, equation (|13|) may not be an accu- 
rate estimate of the self-shielding density threshold at low 
redshifts. Our tests show that simulations that use equa- 
tion (| 14p reproduce the f(Nm,z) accurately to within 10% 
for z > 1 where photoionization is dominant (see Appendix 

ED . 

Although using the relation between the median pho- 
toionization rate and the gas density is a computationally 
efficient way of calculating equilibrium neutral fractions in 
big simulations, it comes at the expense of the informa- 
tion encoded in the scatter around the median photoioniza- 
tion rate at a given density. However, our experiments show 
that the error in f(Nm,z) that results from neglecting the 
scatter in the photoionization rate profile is negligible for 
Nm •> 10 18 cm -3 and less than < 0.1 dex at lower column 
densities (see Appendix lAll) . 

3.4 The roles of diffuse recombination radiation 
and collisional ionization at z = 3 

To study the interplay between different ionizing processes 
and their effects on the distribution of HI, we compare their 
ionization rates at different densities. We start the analysis 
by presenting the results at z = 3 and extend it to other 
redshifts in 333] 



The total photoionization rate profiles shown in the 
right panel of Figure [3] are almost flat at low densities 
and decrease with increasing density, starting at densities 
n H ~ 10~ 4 cm~ 3 . Just below n H = 10~ 2 cm -3 self-shielding 
causes a sharp drop, but the fall-off becomes shallower for 
n H > 10 -2 cm -3 and the photoionization rate starts to in- 
crease at n H > 10 cm" 3 . As shown in Figure[4] the shallower 
fall-off in the total photoionization rate with increasing den- 
sity is caused by RR. The increase in the photoionization 
rate with density at the highest densities on the other hand, 
is an artifact of the imposed temperature for ISM particles 
(i.e., T = 10 4 K) which produces a rising collisional ioniza- 
tion rate with increasing density. As the comparison between 
the UVB and RR photoionization profiles shows (see Figure 
2]), RR om y starts to dominate the total photoionization 
rate at n H > 10~ 2 cm" 3 , where the UVB photoionization 
rate has dropped by more than one order of magnitude and 
the gas is no longer highly ionized. RR reduces the total 
HI content of high-density gas by ~ 20%. Although ion- 
ization rates remain non- negligible at higher densities, they 
cannot keep the hydrogen highly ionized. For instance at 
n H ~ 1cm -3 , a photoionization rate of T ~ 10~ 14 s _1 can 
only ionize the gas by < 20%. 



The shape of the photoionization rate profile produced 
by diffuse RR can be understood by noting that the produc- 
tion rate of RR increases with the density of ionized gas. At 
number densities n H < 10 _2 cm -3 , where the gas is highly 
ionized, the photoionization rate due to recombination pho- 
tons is proportional to the density (i.e., Frr oc n H ). At 
higher densities on the other hand, the gas becomes neu- 
tral. As a result, the density of ionized gas decreases with 
increasing density and the production rate of recombination 
photons decreases. Therefore, there is a peak in the pho- 
toionization rate due to RR around the self-shielding den- 
sity. At very low densities, the superposition of recombina- 
tion photons which have escaped from higher densities be- 
comes dominant and the net photoionization rate of recom- 
bination photons flattens. Note that our simulations may 
underestimate this asymptotic rate because our simulation 
volumes are small compared to the mean free path for ion- 
izing radiation (which is ~ 100 Mpc at z ~ 3). On the 
other hand, the neglect of cosmological redshifting for RR 
will result in overestimation of its photoionization rate on 
large scales. Recombination photons also leak from lower 
densities to self-shielded regions, smoothing the transition 
between highly ionized and highly neutral gas. At high den- 
sities, in the absence of the UVB ionizing photons, RR and 
collisional ionization can boost each other by providing more 
free electrons and ions. 



In Figure [5] we compare hydrogen neutral fraction and 
photoionization rate profiles for different assumptions about 
RR. The hydrogen neutral fraction profile based on a precise 
RT calculation of RR is close to the Case A result at low 
densities (n H < 10 -3 cm -3 ) but converges to the Case B 
result at high densities (n H > 10 _1 cm -3 ). This suggests 
that the neutral fraction profile, though not the ionization 
rate, can be modeled by switching fro m Case A to Case 
B recombination at n„ ~ «H,ssh (e.g.. lAItav et al.ll201ll : 
iMcQuinn et~al1l201lh . 
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Figure 5. The hydrogen neutral fractions {left) and the UVB photoionization rate profiles (right) as a function of density in RT 
simulations with different models for recombination radiation for the L06N256 simulation at z = 3. The red dashed curve shows 
the reference simulation where recombination radiation is modeled self-consistently. The blue solid and green dot-dashed curves show 
simulations in which recombination radiation is substituted by the use of Case A and Case B recombination rates, respectively. The 
curves show the medians and the shaded areas around the medians indicate the 15% — 85% percentiles. HI column densities corresponding 
to each density are shown along the top x-axes. The effect of recombination radiation on the hydrogen neutral fractions is similar to the 
use of Case A recombination at low densities (i.e., n H < 10 -3 cm - 3 ) and to the use of Case B recombination at higher densities (i.e., 
n H > 10 —1 cm -3 ). However, recombination radiation can penetrate into the self-shielded regions, an effect that is not captured by the 
use of Case B recombination. 



3.5 Evolution 

The general trends in the profile of the photoionization rates 
with density and their influences on the distribution of HI 
are not very sensitive to redshift. However, as shown in table 
[2] the intensity and hardness of the UVB radiation change 
with redshift which, in turn, changes the self-shielding den- 
sity. Moreover, as the Universe expands, the average den- 
sity of absorbers decreases and their distributions evolve. 
The larger structures that form at lower redshifts drasti- 
cally change the temperature structure of the gas at low 
and intermediate densities where collisional ionization be- 
comes the dominant process. In the top-left panel of Figure 
[6] the evolution of the hydrogen neutral fraction is illus- 
trated for the L50N512-W3 simulation. As discussed in 33.31 
and Appendix lAll since the photoionization rate profiles are 
converged with box size and resolution, we apply the profiles 
derived from a RT simulation of a smaller box, or a subset 
of the big box at lower redshiftijf], to calculate the neutral 
fractions in this big box. Figure HJ] shows that the neutral 
fraction profiles are similar in shape at high redshifts but 
that at z ^ 1 the profiles are largely different, particularly 
at low hydrogen number densities, due to the evolving col- 
lisional ionization rates. 

The evolution of the collisional ionization rate profiles is 
shown in the bottom-left panel of Figure [6] At z 2 and for 
n H < 10 _2 cm -3 , the collisional ionization rate is not high 
enough to compete with the UVB photoionization rate. At 



8 Strong collisional ionization at low redshifts can change the self- 
shielding. Therefore, for our RT simulation at low redshifts (i.e., 
z < 1) we used a representative sub-volume of the L50N512-W3 
simulation sufficiently large for the collisional ionization rates to 
be converged. 



lower redshifts and for number densities n H < 3x 1CP 3 cm -3 , 
on the other hand, collisional ionization dominate^]. Indeed, 
the median collisional ionization rates are more than 100 
times higher than the UVB photoionization rate at den- 
sities around the expected self-shielding thresholds. Colli- 
sional ionization therefore helps the UVB ionizing photons 
to penetrate to higher densities without being significantly 
absorbed. As a result, self-shielding starts at densities higher 
than expected from equation (I13|l . The signature of colli- 
sional ionization on the hydrogen neutral fraction is more 
dramatic at low densities and partly compensates for the 
lower UVB intensity at z = 0. This results in a flattening of 
f{Nm, z) at column densities Nbj < 10 16 cm -2 as shown in 
the left panel of Figure [1] 

As mentioned above, at low redshifts (e.g., z — 0) the 
collisional ionization rate peaks at densities higher than the 
expected self-shielding threshold against the UVB. As a re- 
sult, the total photoionization rate falls off rapidly together 
with the drop in the collisional ionization rate. Therefore, 
the drop in the hydrogen ionized fraction, and hence the 
resulting free electron density, is much sharper at lower red- 
shifts. This causes a steeper high-density fall-off in the col- 
lisional ionization rate as shown in the bottom-left panel of 
Figure [6] 

The differences between the total photoionization rates 
at different redshifts shown in the top-right panel of Figure 
[6] are caused by the evolution of the UVB intensity and its 
hardness, which affects the self- shielding density thresholds 
(see equation I13[) . On the other hand, as we showed in the 



9 We note that at redshift z < 1 box sizes Lbox ^ 25 comoving 
h~ Mpc are required to fully capture the large-scale accretion 
shocks and to produce converged collisional ionization rates. 
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Figure 6. Evolution of the hydrogen neutral fraction profile {top-left) and various ionization rates as a function of density. Top-right, 
bottom-left and bottom-right panels show, respectively, the total photoionization rates, collisional ionization rates and RR photoionization 
rates. All the HI fractions and collisional ionization rates which are sensitive to both collisional ionization and photoionization are taken 
from the L50N512-W3 simulation. Photoionization rates at z ^ 2 are based on the L06N128 simulation. At lower redshifts (i.e., z = 
and 1), where the box size become important because of the collisional ionization and its effect on changing the self-shielding, we used 
a representative sub-volume of the L50N512-W3 simulation to calculate the photoionization rate profile, with density While the overall 
shape of the UVB photoionization rate profile is similar at different redshifts, the collisional ionization becomes increasingly stronger at 
lower redshifts and strongly reduces the hydrogen neutral fractions at densities n H < 10 -3 cm -3 . 



previous section, the peak of the photoionization rate pro- 
duced by RR tracks the self-shielding density. As a result 
the peak of the RR photoionization rate also changes with 
redshift as illustrated in the bottom-right panel of Figure 
E). 

The filled circles in Figure[7|indicate the UVB photoion- 
ization rate versus the number density at which the RR pho- 
toionization rate peaks. The self-shielding density expected 
from the Jeans scaling argument (equation [TH) is also shown 
(green dotted line). The peaks in the RR photoionization 
rate in RT simulations follow this expected scaling for z ^ 1. 
However, the z = result deviates from this trend since col- 
lisional ionization affects the self-shielding density threshold, 
a factor that is not captured by equation (|13p. 

As a result of the RR photoionization rate peaking 
around the self-shielding density threshold, the transition 
between highly ionized and nearly neutral gas becomes more 
extended at all redshifts. To illustrate the smoothness of this 
transition, the densities at which the median hydrogen neu- 
tral fractions are 10 -2 and 0.5 are shown in Figure [7] with 
open circles and squares, respectively. The densities at which 



the hydrogen neutral fraction is 10 -2 are slightly higher than 
the densities at which the RR photoionization rate peaks 
(filled circles). The evolution agrees with the trend expected 
from the Jeans scaling argument and the self-shielding den- 
sity. The exception is again z — 0, where the large collisional 
ionization rate at number densities n H ~ 10 -3 — 10 -2 cm -3 
shifts the transition to neutral fraction of 0.5 to densities 
that are « 1 dex higher. However, the relation between the 
photoionization rate and the density still follows the slope 
expected from equation (|13l) . 



It is interesting to note that the UVB spectral shape 
at z = 4 is slightly harder than at z = 1 while the UVB 
intensities at these two redshifts are similar. This results in 
a deeper penetration of ionizing photons at z = 4. Conse- 
quently, the densities corresponding to the indicated neutral 
fractions (i.e., 10 -2 and 0.5) at z = 1 are lower than their 
counterparts at z = 4. 
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Figure 7. The photoionization from recombination radiation 
peaks at the expected self-shielding density and smooths the tran- 
sition between highly ionized gas and self-shielded gas. This is il- 
lustrated by showing different characteristic densities for various 
redshifts. Filled circles show the density at which recombination 
radiation peaks, while open circles and squares show the HI num- 
ber density corresponding to hydrogen neutral fractions of 10 — 2 
and 0.5, respectively. The self-shielding density threshold for a 
given photoionization rate expected from the Jeans scaling ar- 
gument (equation I13H is also indicated by the green dotted line. 
The Jeans scaling argument works well, except at z = when 
collisional ionization is important. 



4 CONCLUSIONS 

We combined a set of cosmological hydrodynamical simula- 
tions with an accurate RT simulation of the UVB radiation 
to compute the HI column density distribution function and 
its evolution. We ignored the effect of local sources of ioniz- 
ing radiation, but we did include a self-consistent treatment 
of recombination radiation. 

Our RT results for the distribution of photoionization 
rates at different densities are converged with respect to the 
simulation box size and resolution. Therefore, the resulting 
photoionization rate can be expressed as a function of the 
hydrogen density and the UVB. We provided a fit for the 
median total photoionization rate as a function of density 
that can be used with any desired UVB model to take into 
account the effect of HI self-shielding in cosmological simu- 
lations without the need to perform RT. 

The CDDF, /(Vhi, z), predicted by our RT simulations 
is in excellent agreement with observational constraints at 
all redshifts (z = — 5) and reproduces the slopes of the ob- 
served /(JVhi, z) function for a wide range of HI column den- 
sities. At low HI column densities, the CDDF is a steep func- 
tion which decreases with increasing Njjx before it flattens at 
Nm > 10 18 cm" 2 due to self-shielding. At AT H i > 10 21 cm" 3 
on the other hand, f(Nm,z) is determined mainly by the 
intrinsic distribution of total hydrogen and the H2 fraction. 

We showed that the Nm — n H relationship can be ex- 
plained by a simple Jeans scaling. This argument assumes HI 
absorbers to be sel f-gravitating s ystems close to local hydro- 
static equilibrium (Schavc 2001al) and to be either neutral or 
in photoionization equilibrium in the presence of an ionizing 
radiation field. However, at z — the analytic treatment un- 



derestimates the self-shielding density threshold due to its 
neglect of collisional ionization. 

The high HI column density end of the predicted 
/(Vhi, z) evolves only weakly from z — 5 to z = 0, consistent 
with observations. In the Lyman limit range of the distri- 
bution function, the slope of /(Vhi, z) remains the same at 
all redshifts. However, at z > 3 the number of absorbers in- 
creases with redshift as the Universe becomes denser while 
the UVB intensity remains similar. At lower redshifts, on the 
other hand, the combination of a decreasing UVB intensity 
and the expansion of the Universe results in a non-evolving 
/(Vhi, z). In contrast, the number of absorbers with lower 
HI column densities (i.e., the Lya forest) decreases signifi- 
cantly from z ~ 3. We showed that this results in part from 
the stronger collisional ionization at redshifts z < 1, which 
compensates for the lower intensity of the UVB. The increas- 
ing importance of collisional ionization is due to the rise in 
the fraction of hot gas due to shock-heating associated with 
the formation of structure. 

The inclusion of diffuse recombination radiation 
smooths the transition between optically thin and thick 
gas. Consequently, the transition to highly neutral gas is 
not as sharp as what has been assumed in some previ- 
ous works (e . g., Nagamine et all 120101 ; lYaiima et alj l201ll : 
iGoerdt et alJl2012i ). For instance, the difference in the gas 
density at which hydrogen is highly ionized (i.e., n H1 /n H < 
0.01) and the density at which gas is highly neutral (i.e., 
n m/ n H ?i 0-5) is more than one order of magnitude (see 
Figure [7]). As a result, assuming a sharp self-shielding den- 
sity threshold at the density for which the optical depth of 
ionizing photons is ~ 1, overestimates the resulting neutral 
hydrogen mass by a factor of a few. 

Our simulations adopted some commonly used approxi- 
mations (e.g., neglecting helium RT effects, using a gray ap- 
proximation in order to mimic the UVB spectra, neglecting 
absorption by dust and local sources of ionizing radiation). 
Our tests show that most of those approximations have neg- 
ligible effects on our results. But there are some assumptions 
which require further investigation. For instance, the pres- 
ence of young stars in high-density regions could change the 
HI CDDF, especially at high HI column densities through 
feedback and emission of ionizing photons. Indeed, we will 
show in Rahmati et al. (in prep.) that for very high column 
densities the ionizing radiation from young stars can reduce 
the /(V m , z) by 0.5-1 dex. 
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Table Al. The best-fit parameters for equation IIAH at different 
redshifts based on RT results in the L06N128 simulation. 



Redshift 


log [n () ] (cm 3 ) 


ai 


"2 


P 


1-/ 


z = 1-5 


log [™H,SSh] 


-2.28 


-0.84 


1.64 


0.98 


2 = 


-2.94 


-3.98 


-1.09 


1.29 


0.99 


2 = 1 


-2.29 


-2.94 


-0.90 


1.21 


0.97 


2 = 2 


-2.06 


-2.22 


-1.09 


1.75 


0.97 


2 = 3 


-2.13 


-1.99 


-0.88 


1.72 


0.96 


2 = 4 


-2.23 


-2.05 


-0.75 


1.93 


0.98 


2 = 5 


-2.35 


-2.63 


-0.57 


1.77 


0.99 



APPENDIX A: PHOTOIONIZATION RATE AS 
A FUNCTION OF DENSITY 

Al Replacing the RT simulations with a fitting 
function 



In H3.3I we demonstrated that the median of the simulated 
relation between the total photoionization rate, Tphot, and 
density is converged with respect to resolution and box size. 
We used this result and provided fits to the median of this 
relation. We have exploited these fits to compute the neutral 
hydrogen fraction in cosmological simulations under the as- 
sumption of ionization equilibrium (see Appendix lA2[) . with- 
out performing the computationally demanding RT. In this 
section, we discuss the accuracy of these fits. 

The left panel of Figure Q] shows that using the median 
photoionization rates produces an HI CDDF in very good 
agreement with the HI CDDF obtained from the correspond- 
ing RT simulation (orange solid curve) at Nm > 10 18 cm -2 . 
However, there is a small systematic difference at lower col- 
umn densities. One may think that this small difference is 
caused by the loss of information contained in the scatter 
in the photoionization rates at fixed density. We tested this 
hypothesis by including a log-normal random scatter around 
the median photoionization rate consistent with the scatter 
exhibited by the RT result. However, after accounting for 
the random scatter, the f(Nm,z) is slightly overproduced 
compared to the full RT result at nearly all HI column den- 
sities. 

We exploit the insensitivity of the shape of the Tphot- 
density relation to the redshift, and propose the following 
fit to the photoionization rate, Tphot, 



Fuvb 



(1-/) 



1 + 



Us. 

no 



+ .f 



1 + 



no 



(Al) 



where Fuvb is the photoionization rate due to the ioniz- 
ing background, and no, Qi, as, and j3 are parameters of 
the fit. The best-fit values of these parameters are listed 
in Table IA1I and the photoionization rate-density relations 
they produce are compared with the RT simulations at red- 
shifts 2 = and z = 4 in Figure IA1I At all redshifts, the 
best-fit value of no is almost identical to the self-shielding 
density threshold, nH,ssh, defined in equation 1)130 . and the 
characteristic slopes of the photoionization rate-density re- 
lation are similar. This suggest that one can find a sin- 
gle set of best-fit values to reproduce the RT results at 
z > 1. The corresponding best-fit parameter values are (see 
also equation [14)) a x = -2.28 ± 0.31, a 2 = -0.84 ± 0.11, 



n = (1.003 ± 0.005) x n H ,ssh, P = 1-64 ± 0.19 and / = 
0.02 ± 0.0089. 

In the right panel of Figure [1] the ratio between the 
HI CDDF calculated using the fitting function presented in 
equation l)14fl and the RT based HI CDDF (i.e., calculated 
using the median of the photoionization rate-density relation 
in the RT simulations) is shown for the L50N512-W3 and 
at z — 0, 2 and 4. This illustrates that the fitting function 
reproduces the RT results accurately, except at z — 0. As 
explained in [|3] this is expected since at low redshifts colli- 
sional ionization affects the self-shielding and the resulting 
photoionization rate-density profile. However, a separate fit 
can be obtained using converged RT results at z = 0. 



A2 The equilibrium hydrogen neutral fraction 

In this section we explain how to derive the neutral fraction 
in ionization equilibrium. Equating the total number of ion- 
izations per unit time per unit volume with the total number 
of recombinations per unit time per unit volume, we obtain 



nm Ttot = a a n c nail, 



(A2) 



where n HI , n e and the number densities of neu- 

tral hydrogen atoms, free electrons and protons, respectively. 
Ttot is the total ionization rate per neutral hydrogen atom 
and cya is the Case A recom bination rat for wh ich we use 
the fitting function given bv lHui fc Gnedinl (| 19971 ): 



a A = 1.269x10" 



(1 + (A/0.522) 047 ) 



0.47-v 1-923 



cm 3 s"\ (A3) 



where A = 315614/T. 

Defining the hydrogen neutral fraction as the ratio be- 
tween the number densities of neutral hydrogen and total 
hydrogen, n = n HI /n H , and ignoring helium (which is an 
excellent approximation, see Appendix ID2)) . we can rewrite 
equation (|A3[) as: 



V Ttot = a A (1 - v) n u- 



(A4) 



Furthermore, we can assume that the total ionization rate, 
Ttot, consists of two components: the total photoionization 



10 The use of Case B is more appropriate for n H > rtjj sSh- How- 
ever, we assume the photoionization due to RR is included in 
Ttot, e.g., by using the best-fit function that is presented in 
equation [14] Therefore, Case A recombination should be adopted 
even at high densities. 
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Figure 1. Left: The ratio between the HI CDDF calculated using the RT based Tp^ot ~ " H relationship and the actual RT results for 
the L06N128 simulation in the presence of the UVB and diffuse recombination radiation, at Z = 3. The orange solid line shows the 
result of using the median photoionization rate-density profile predicted by the RT simulation. The blue dashed curve shows the result 
of including the scatter around the median in the calculations. Right: HI CDDFs calculated using the Fpi lot — ra H fitting function (i.e., 
equation I14H are compared to the HI CDDFs for which the actual rp hot -density relation from the RT simulations are used. Blue and 
green curves are for 2 = and z = 2 respectively and the red curve is for 2 = 4. All the CDDFs are for the L50N512- W3 simulation and 
in the presence of the UVB and diffuse recombination radiation. 
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Figure Al. Comparisons between the total photoionization rates as a function of density based on the RT simulations and best-fit 
functions at z = 4 and 2 = are shown in the left and right panels, respectively. In each panel, the RT result is shown with the orange 
solid curve. The best fit to the RT result at a given redshift (equation lAll and Table |AH is shown with the blue dashed curve and the 
best fit to the RT results at z = 1 — 5 (equation 1 14H is shown with the purple dotted curve. As shown in the right panel, because of the 
impact of collisional ionization on self-shielding, the low redshift photoionization curve (the blue dashed curve) deviates from the best 
fit to the results at higher redshifts (the purple dotted curve). 



rate, Tphot, and the collisional ionization rate, Fcoi: 

Ttot = Fphot + Tcoi, (A5) 

where Fcoi = At (1 — ??) n H . The photoionization rate can be 
expressed as a function of density using equation (|14[) . For 
At, whi ch depends on l y on temperature, we use a relation 
given in iTheuns et aTJ l| 1998ft : 



At = 1.17 x ID" 10 ^exp(-157809/r) ^-t } 



i + v/ryicF 



with A — ola + A T , B = 2a A + 
which gives: 

Vb^ 



+ At and C — qa 



B 



V - 



4AC 



2A 



(A8) 



We can now rearrange equation (|A4I> as a quadratic 
equation: 



A<n* -B ri+C = 0, 



(A7) 



Using the last equation one can calculate the equilib- 
rium hydrogen neutral fraction for a given n H and temper- 
ature. 
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APPENDIX B: THE EFFECTS OF BOX SIZE, 
COSMOLOGICAL PARAMETERS AND 
RESOLUTION ON THE HI CDDF 

The size of the simulation box may limit the abundance and 
the density of the densest systems captured by the simula- 
tion. In other words, very massive structures, which may be 
associated with the highest HI column densities, cannot be 
formed in a small cosmological box. Indeed, as shown in the 
top panels of Figure |BT1 one needs to use cosmological boxes 
larger than > 25 comoving /i -1 Mpc in or der to achieve con - 
vergence in the HI distribution (see also lAltav et al.ll201lf ). 
On the other hand, the bottom-right panel of Figure IB1I 
shows that changing the resolution of the cosmological sim- 
ulations also affects f(Nm,z), although the effect is small. 

The adopted cosmological parameters also affect the gas 
distribution and hence the HI CDDF. For instance, one ex- 
pects that the number of absorbers at a given density varies 
with the density parameter f2b, and the root mean square 
amplitude of density fluctuations as- The bottom- left panel 
of Figure IB 1 1 shows the ratio of column densities in simula- 
tions assuming WMAP 7-year and 3-year parameters. The 
ratio is only weakly dependent on the box size of the simu- 
lation and its resolution. This motivates us to use this ratio 
to convert the HI CDDF between the two cosmologies for 
all box sizes and resolutions (at any given redshift). While 
this is an approximate way of correcting for the difference 
in the cosmological parameters, it does not affect the main 
conclusions presented in this work (e.g., the lack of evolution 
of f(N m ,z)). 



APPENDIX C: RT CONVERGENCE TESTS 
CI Angular resolution 

The left-hand panel of Figure [CTI shows the dependence of 
photoionization rates on the adopted angular resolution, i.e., 
the opening angle of the transmission cones 4tt/Ntc- The 
photoionization rates are converged for iVxc = 64 (our fidu- 
cial value) or higher. 

C2 The number of ViP neighbors 

The right panel of Figure I C 1 1 shows the dependence of the 
photoionization rates on the number of SPH neighbors of 
ViPs. As discussed in H2.21 ViPs distribute the ionizing pho- 
tons they absorb among their NGBviP nearest SPH neigh- 
bors. The larger the number of neighbors, the larger the 
volume over which photons are distributed, and the more 
extended is the transition between highly ionized and self- 
shielded gas. The photoionization rates converge for < 5 ViP 
neighbors (our fiducial value is 5). 

C3 Direct comparison with another RT method 

lAltav et al.l (|201 lh used cosmological s imulations from th e 
reference model of the OWLS project (jSchave et al.l 12010^) . 
i.e., a simulation run with the same hydro code as we used 
in this work, to investigate the effect of the UVB on the 
HI CDDF at z — 3. However, they employed a ray-tracing 
method very different from the RT method we use here. 



Furthermore, they did not explicitly treat the transfer of re- 
combination radiation. In Figure [C2l we compare one of our 
UVB photoioniza tion rate profile^ 1 1 ! with the photoioniza- 
tion rate found bv lAltav et al.l |201ll ) in a similar simulation. 
The overall agreement is very good, but the comparison also 
reve als important diffe rences. 

lAltav et al.l (|201 lh calculate the average optical depth 
around every SPH particle within a distance of 100 proper 
kpc, assuming the UVB is unattenuated at larger distances. 
Then, they use this optical depth to calculate the attenua- 
tion of the UVB photoionization rate for every particle. This 
procedure may underestimate the small but non-negligible 
absorption of UVB ionizing photons on large scales. Indeed, 
by tracing the self-consistent propagation of photons inside 
the simulation box, we have found that the UVB photoion- 
ization rate decreases gradually with increasing density up 
to the density of self-shielding. However, we note that the 
small differences betwe en our UVB phot oionization rates 
and those calculated bv lAltav et all l|201ll) at densities be- 
low the self-shielding, become slightly smaller by increasing 
the angular resolution in our RT calculations (see the left 
panel of Figure ICTjl . 

APPENDIX D: APPROXIMATED PROCESSES 
Dl Multifrequency effects 

As discussed in £12.31 in our RT simulations we have treated 
the multifrequency nature of the UVB radiation in the gray 
approximation (see equation [4]). This approach does not 
capture the spectral hardening which is a consequence of 
variation of the absorption cross-sections with frequencies. 
We tested the impact of spectral hardening on the HI frac- 
tions by repeating the L06N128 simulation at z — 3 with 
the UVB using 3 frequency bins. We used energy intervals 
[13.6 - 16.6], [16.6 - 24.6] and [24.6 - 54.4] eV and assumed 
that photons with higher frequencies are absorbed by He. 
The result is illustrated in the top section of the left panel 
in Figure ID11 by plotting the ratio between the resulting 
hydrogen neutral fraction, r/, and the same quantity in the 
original simulation that uses the gray approximation. This 
comparison shows that the simulation that uses multifre- 
quency predicts hydrogen neutral fractions < 10% lower at 
low densities (i.e., n H < 10 _4 cm~ 3 ). This does not change 
the resulting /(JVhi, z) noticeably at the column densities of 
interest here. 

The spectral hardening captured in the simulation with 
3 frequency bins is illustrated in the right panel of Figure 
ID1I This figure shows the fractional contribution of differ- 
ent frequencies to the total UVB photoionization rate as a 
function of density. The red solid curve shows the contribu- 
tion of the bin with the lowest frequency and drops at the 
self-shielding density threshold. On the other hand, the frac- 
tional contribution of the hardest frequency bin increases at 
higher densities, as shown with the blue dashed curve. De- 
spite the differences in the fractional contributions to the to- 
tal UVB photoionization rate, the absolute photoionization 



Note that in our simulations the UVB photoionization rate is 
converged with the box size and the resolution as shown in £|3. 31 
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Figure Bl. The relative changes in the HI CDDF using different resolutions, box sizes and cosmologies in the presence of the UVB 
and diffuse recombination radiation. The top-left panel shows the effect of box size on /(JVhIi z) for a fixed resolution at z = 3, where 
the orange solid (blue dashed) curve shows the difference between using a box size of L = 25 (50) comoving h — 1 Mpc and a box size 
of L = 100 comoving ft - *Mpc. The top-right panel shows the same effect but for smaller box sizes: the orange solid (blue dashed) 
curve shows the difference between using a box size of L = 6 (12) comoving ft - 1 Mpc and a box size of L = 25 comoving ft - 1 Mpc. 
The bottom-left shows the effect of using a cosmology consistent with WMAP 3-year results instead of using a cosmology based on the 
WMAP 7-year constraints. The orange solid and blue dashed curves show this effect for simulations with box sizes of L = 6 and 25 
comoving ?t — 1 Mpc, respectively. The bottom-right panel shows the effect of resolution. 
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Figure CI. The UVB photoionization rate is converged for our adopted angular resolution, i.e., Ntc = 64, as shown in the left panel 
and our adopted number of ViP neighbors, i.e., NGBvip = 5, as shown in the right panel. Photoionization rate profiles are shown for the 
L06N128 simulation in the presence of the UVB radiation where the Case A recombination is adopted. The curves show the medians 
and the shaded areas around them indicate the 15% — 85% percentiles. 
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Figure C2. Left: Median UVB photoionization rate as a fun ction of dens i ty at z = 3 using different RT methods. The red dashed curve 
shows the results based on the method that has been used in lAItav et al . (2011) and the blue solid curve shows the result of this work. 
Right: The HI CDDF of the L25N512 si mulation at 2 = 3 using different RT methods and without RR. The red dashed curve shows 
the ratio between the H I CDDF given in | Altav et al] (|201 lh and our results. This comparison shows that despite the overall agreement 
between our results and lAltav et al.l ll201lh . there are some important differences. 
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Figure Dl. Spectral hardening and multifrequency treatment do not change the HI distribution significantly. Left: The ratio between 
hydrogen neutral fractions, r), obtained by using 3 frequency bins and by using the gray approximation is shown in the top-left apnel. 
The ratio between hydrogen neutral fractions resulting from a simulation with 4 frequency bins and explicit He treatment and the same 
quantity using the gray approximation and without explicit He treatment is shown in the bottom-left panel. The vertical lines with 
different lengths indicate the 15% — 85% percentiles. Right: The fractional contribution of different frequency bins to the total UVB 
photoionization rate for the simulations with 3 frequency bins. All the RT calculations are performed using the L06N128 simulation at 
z = 3 in the presence of the UVB and assuming Case A recombination. 



rates drop rapidly at densities higher that the self-shielding 
threshold for all frequency bins. 



D2 Helium treatment 

A simplifying assumption frequently used in RT simulations 
which aim to calculate the distribution of neutral hydro- 
gen is to ignore helium in t h e ionization processe s (e.) 



Faucher-Giguere et alj 120091 ; iMcQuinn fc SwitzeJ [201 



Altav et al. 201 lh . We adopted the same assumption in our 



RT calculations which implies that we implicitly assumed 
the ionization state of neutral helium and its interaction 
with free electrons to be similar to the trends followed by 
neutral h ydrogen. This has been shown to be a good as- 
sumption (Ostcrbrock & Fcrland 2 0061 : IMcQuinn fc Switzen 
120101 : iFriedrich et alll2012t ). Nevertheless, we tested the va- 
lidity of our approximate helium treatment by repeating the 



L06N128 simulation at z = 3 with the UVB using 4 fre- 
quency bins and an explicit He treatment. The first three 
frequency bins are identical to the bins used in the previous 
section (i.e., [13.6 - 16.6], [16.6 - 24.6] and [24.6 - 54.4] eV) 
and the last bin is chosen to cover higher frequencies which 
are capable of Hell ionization. We adopted a helium mass 
fraction of 25% and a Case A recombination rate. The ra- 
tio between the resulting hydrogen neutral fraction and the 
same quantity when a single frequency is used and helium is 
not treated explicitly is illustrated in the bottom-left panel 
of Figure IDTI The hydrogen neutral fractions are very close 
in the two simulations. However, the simulation with mul- 
tifrequency and explicit He treatment results in hydrogen 
neutral fractions that are < 10% higher at low densities 
(i.e., n H < 10~ 4 cm -3 ). This difference is barely noticeable 
in the comparison between the two HI CDDFs (not shown). 
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